import numpy as np
from numpy import pi, sin, exp, log
import plotly.graph_objects as go
import matplotlib.pyplot as plot
def f(x,y):
return 10*((x)**5 + pi)*sin(log(((y)**4)+1))/exp(0.8*x**2) + x*5 + 1*(1 - x)**2 - y*2
global_epsilon = 0.000000001 # argument increment for derivative
radius = 8 # working plane radius
centre = (global_epsilon, global_epsilon) # centre of the working circle
arr_shape = 100 # number of points processed / 360
step = radius / arr_shape # step between two points
def rotate_vector(length, a):
return length * np.cos(a), length * np.sin(a)
def derivative_y(epsilon, arg):
return (f(arg, epsilon + global_epsilon) - f(arg, epsilon)) / global_epsilon
def derivative_x(epsilon, arg):
return (f(global_epsilon + epsilon, arg) - f(epsilon, arg)) / global_epsilon
def calculate_flip_points():
flip_points = np.array([0, 0])
points = np.zeros((360, arr_shape), dtype=bool)
cx, cy = centre
for i in range(arr_shape):
for alpha in range(360):
x, y = rotate_vector(step, alpha)
x = x * i + cx
y = y * i + cy
points[alpha][i] = derivative_x(x, y) + derivative_y(y, x) > 0
if not points[alpha][i - 1] and points[alpha][i]:
flip_points = np.vstack((flip_points, np.array([alpha, i - 1])))
return flip_points
def pick_estimates(positions):
vx, vy = rotate_vector(step, positions[1][0])
cx, cy = centre
best_x, best_y = cx + vx * positions[1][1], cy + vy * positions[1][1]
for index in range(2, len(positions)):
vx, vy = rotate_vector(step, positions[index][0])
x, y = cx + vx * positions[index][1], cy + vy * positions[index][1]
if f(best_x, best_y) > f(x, y):
best_x = x
best_y = y
for index in range(360):
vx, vy = rotate_vector(step, index)
x, y = cx + vx * (arr_shape - 1), cy + vy * (arr_shape - 1)
if f(best_x, best_y) > f(x, y):
best_x = x
best_y = y
return best_x, best_y
def gradient_descent(best_estimates, is_x):
derivative = derivative_x if is_x else derivative_y
best_x, best_y = best_estimates
descent_step = step
value = derivative(best_y, best_x)
while abs(value) > global_epsilon:
descent_step *= 0.95
best_y = best_y - descent_step \
if derivative(best_y, best_x) > 0 else best_y + descent_step
value = derivative(best_y, best_x)
return best_y, best_x
def find_minimum():
return gradient_descent(gradient_descent(pick_estimates(calculate_flip_points()), False), True)
def draw_chart(point):
n_points = 400
x_range = (-5, 5)
y_range = (-5, 5)
X = np.linspace(x_range[0], x_range[1], n_points)
Y = np.linspace(y_range[0], y_range[1], n_points)
Z = np.zeros((n_points, n_points))
for i, x in enumerate(X):
for j, y in enumerate(Y):
Z[j][i] = f(x, y)
x, y, z = minimum
scatter_dict = {
'x':[float(x)],
'y':[float(y)],
'z':[float(z)],
'mode':'markers',
'marker':go.scatter3d.Marker(color='green'),
}
fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y), go.Scatter3d(**scatter_dict)])
fig.update_layout(title='f(x,y) 3D plot', autosize=False,
width=1000, height=1000)
fig.show()
min_x, min_y = find_minimum()
minimum = (min_x, min_y, f(min_x, min_y))
print(minimum)
draw_chart(minimum)
(-0.05761390345066797, 3.284635397725483, -37.027090239286075)